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Abstract 

Boundary scattering of thermal phonons in thin solid films is typically analyzed using Fuchs- 
Sondheimer theory, which provides a simple equation to calculate the reduction of thermal conduc¬ 
tivity as a function of the film thickness. However, this widely-used equation is not applicable to 
highly anisotropic solids like graphite because it assumes the phonon dispersion is isotropic. Here, 
we derive a generalization of the Fuchs-Sondheimer equation for solids with arbitrary dispersion 
relations and examine its predictions for graphite. We find that the isotropic equation vastly over¬ 
estimates the boundary scattering that occurs in thin graphite films due to the highly anisotropic 
group velocity, and that graphite can maintain its high in-plane thermal conductivity even in thin 
films with thicknesses as small as ten nanometers. 
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Thermal transport in thin solid hlms with thicknesses from tens of nanometers to mi¬ 
crometers is a topic of considerable importance/^^ with snch hlms being used in applications 
ranging from quantum well lasers to LED lightingT^^ In these thin hlms, the phonon mean 
free paths (MFPs) are comparable to the him thickness, resulting in boundary scattering that 
reduces the thermal conductivity and inhibits heat dissipation. This thermal management 
problem is presently an important challenge in many devices such as GaN transistors.®’®’® 

Heat transport along the in-plane direction of thin hlms is typically described using the 
Fuchs-Sondheimer equation, which is an analytical solution of the Boltzmann transport 
equation (BTE) for thin hlms with partially specular and partially dihuse walls. This 
equation was originally derived for electron transport^®’^^ and later extended to phonon 
thermal transport assuming an average phonon MFP, enabling the calculation of thermal 
conductivity as a function of the him thickness.Mazumder and Majumdar used a Monte- 
Carlo method to study the phonon transport along a silicon thin him including phonon 
dispersion and polarizations.^^ 

While the Fuchs-Sondheimer equation is in wide use, it cannot be applied to anisotropic 
transport in its typical form because it assumes the crystal of interest is isotropic. However, 
there are many situations in which boundary scattering occurs in thin anisotropic hlms, 
with the most familiar example being thin graphite hlms. Such hlms have been studied 
experimentally,^® and few-layer graphene hlms have been investigated as in heat spreaders 
for GaN transistors.® Provided that the hlms are sufficiently thick that phonon dispersion 
modihcations in the cross-plane direction can be neglected, mathematically describing ther¬ 
mal transport in anisotropic thin hlms requires a Fuchs-Sondheimer equation that is valid for 
any crystal, regardless of its anisotropy. Surprisingly, despite the simplicity of the derivation, 
no such equation has been reported to the best of our knowledge. 

Here, we report the generalization of Fuchs-Sondheimer theory to crystals with arbitrary 
anisotropies. We hnd that highly anisotropic solids with small group velocities along certain 
crystallographic directions experience minimal boundary scattering because the thermal 
conductivity reduction depends only on the component of the MFP normal to the boundary 
rather than the overall MFP. As a result, thin hlms of anisotropic crystals like graphite 
maintain their high thermal conductivity even as the him thickness becomes very small. 
This observation has important implications for heat spreading in electronic devices and the 
thermal conductivity of graphite foams. 
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We begin by considering the steady BTE in a thin film under the relaxation time ap¬ 
proximation. We let the y direction represent the cross-plane direction and assume that a 
thermal gradient exists along the direction. Then, the BTE is given by: 

dk , dk k-fo{T{z,t)) 

= - ( 1 ) 
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where k is the desired distribution function, Vy^ and Vzk are the components of the phonon 
group velocity along the y and directions, Tk is the phonon relaxation time, and k is the 
phonon wavevector in phase space. This equation is solved for the thin film by letting the 
small perturbation = f — fo(T{z,t)), yielding: 


A , _ A dT 

■“^j/k Q “T g]i i^k 


dy 


dT dz 


( 2 ) 


This result assumes that dg\^/dz ~ 0 and that the length scale of the thermal gradient 
is long compared to the phonon MFPs so that a temperature gradient can be defined. 
Equation 2 is a one-dimensional ODE with the boundary conditions gkk < ^lU = d) = 
gkk > 0 , 1 / = 0) = 0 corresponding to thermalizing, diffuse scattering, meaning that the 
phonon distribution emerging from the wall is a thermal distribution at the local boundary 
temperature. The general solution of this equation is well-known and given in Ref. as: 



( 3 ) 
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where g = ///A^^k, k = d/Ay\^ is the nondimensional thickness defined by the y compo¬ 
nent of the MFP, and S'o(; 2 ) = — A^k^^- The thermal conductivity can be obtained by 
substituting this solution into the expression for heat flux: 
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where V is the crystal volume. After evaluating this expression, the thermal conductivity 
of the thin film along the direction is identified as: 
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where the intrinsic thermal conductivity K,zz,]i = C'k'y^k''’k- S{^y), which we term the 
boundary scattering function, is the generalized version of the Fuchs-Sondheimer expression 
describing the reduction in thermal conductivity that occurs due to boundary scattering. 
For reference, the traditional Fuchs-Sondheimer expression obtained by integrating over all 
solid angles in a spherically symmetric Brillouin zone is^^: 

S(?) = 1 - ^ [1 - 4B3(?) + 4B3(0] (7) 

where i?n(0 is the exponential integral function and ^ = d/A^ where is the isotropic 
MFP that depends only on phonon frequency. 

Comparing Eqs. 6 and 7, we see that the traditional Fuchs-Sondheimer expression aver¬ 
ages out the solid angle dependence in phase space so that the hnal expression depends only 
on the overall MFP. The generalized expression does not make any assumptions about the 
symmetry of the crystal, and so the expression depends on the components of the MFPs 
along the crystal axis normal to the boundary. If this MFP component is different than the 
value along the in-plane axis, the actual in-plane thermal conductivity obtained with Eq. 6 
may be very different than the prediction of the traditional expression. 

To examine this prediction, we consider transport along the basal plane of a thin him 
of graphite. We use a phonon dispersion calculated from an optimized Tersoff potential 
provided by Lucas Lindsay.^® The actual phonon-phonon relaxation times of graphite are 
not readily available, and thus we must assume a general form to proceed. In a previous 
work for which we modeled the graphite dispersion with an anisotropic Debye model, 
we found that isotropic relaxation times were able to explain the thermal conductivity 
anisotropy of graphite. In this work for which we use the actual dispersion, we hnd that 
these same relaxation times yield Kyy = 2000 W/mK and Kzz = 230 W/mK. Given that 
the in- and cross-plane thermal conductivities of graphite are 2000 W/mK and 6 W/mK,^° 
respectively, this calculation indicates that these relaxation times cannot explain the large 
thermal anisotropy of graphite. 

We now must hnd the relaxation times that best reproduce the following experimental 
observations. First, the thermal conductivity along each crystal axis is well-known.Second, 
recent computational and experimental reports indicate that phonon MFPs along the c-axis 
are on the order of several hundred nanometers.We satisfy these constraints by assuming 
relaxation times of the form^^: 
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T,{u)-^ = AJ^y ( 8 ) 

\^max,p J 

where Tp{u) is the relaxation time for branch p, Ap is a constant for each branch, u 
is the phonon frequency, and 0Jmax,p is the maximum phonon frequency for each branch. 
Because identifying phonon polarizations off high-symmetry points in the Brillouin zone is 
challenging, branches are determined by sorting the frequencies obtained from the dynamical 
matrix from low to high. 

We can satisfy the constraints by changing Ap for specihc branches even though the 
relaxation times remain isotropic because the different branches contribute differently to 
thermal conductivity along each crystal axis. We take Ap = Apo = 8 x 10^® s“^ for all 
branches except 1, 2, 3, and 5. Branches 1 and 2 primarily carry heat along the c-axis and 
must have smaller relaxation times to yield the correct c-axis thermal conductivity. We 
hnd that Ap^i = Ap^o/lOO and Ap ^2 = ^p,o/20. Branch 3 must have shorter relaxation times, 
yielding Ap^^ = Ap^o/20. We also cap the maximum relaxation times at 250 ps; otherwise the 
c-axis MFPs are too long compared to experiment. Branch 5, the LA branch along the ab 
axis, must have longer relaxation times to reproduce the high ab-axis thermal conductivity, 
Ap^z = 3 X Ap 0 , with a maximum relaxation time of 500 ps. Finally, the minimum relaxation 
time of all branches is enforced to be 5 ps so that the MFPs are not unphysically short. These 
relaxation times yield Kab = 1500 W/mK and Kc = 6.9 W/mK, in reasonable agreement with 
experiment. We emphasize that these relaxation times are not necessarily those of actual 
graphite, but the thermal transport properties obtained using them are sufficiently close to 
the experimental values that they can be used for further analysis. 

The accumulative thermal conductivity versus the component of the MFP along the a, 
b, and c crystal axes is shown in Fig. 1, showing that the chosen relaxation times result in 
phonons along the ab-axis having MFPs on the order of microns to tens of microns, as might 
be expected due to the high basal plane thermal conductivity. On the other hand, phonons 
along the c-axis have MFPs up to 1 micron and typically in the hundreds of nanometers 
range, much shorter than those along the ab axis but still considerably longer than prior 
estimates.The calculated MFPs are, however, of the same order as those recently observed 
in experiments^ and calculated with molecular dynamics. 

The large difference in MFP between the ab and c axes affects thermal transport in thin 
hlms because the strength of boundary scattering only depends on the component of MFP 
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FIG. 1: Normalized, accumulated thermal conductivity versus component of MFP along the trans¬ 
port direction for each crystal axis using the actual graphite phonon dispersion and semi-empirical 
relaxation times. The c-axis MFPs are on the order of hundreds of nanometers, in agreement with 
experiment and considerably shorter than those along the ab-axis.^^ 


normal to the boundary. In the case of transport in a thin him along the ab-axis, the far 
shorter c-axis MFP leads to only minimal boundary scattering even for him thicknesses as 
small as 10 nm. 

The ab-axis thermal conductivity versus him thickness obtained using Eq. 6 is shown in 
Fig. 2a. The ab-axis thermal conductivity remains close to 1000 W/mK even at thicknesses 
of around 10 nm, and nearly reaches its ultimate value at a thickness of around 1 micron. 
These observations are quite unexpected if one considers that the MFPs along the ab-axis 
are tens of microns as in Fig. 1. 

To place these calculations in perspective, we seek to calculate the thermal conductivity 
of an equivalent isotropic crystal with the same ab-axis MFPs. However, we were unable 
to hnd a simple way to convert the highly anisotropic graphite dispersion to an equivalent 
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FIG. 2: (a) Actual ab-axis thermal conductivity using Eq. 6 (solid line) versus film thickness and 
comparison using the same equation with the ab-axis rather than c-axis component of the MFPs in 
the boundary scattering function (dashed line). The actual thermal conductivity remains high even 
for thicknesses as small as 10 nm due to the small MFP component along the thickness direction 
of the film, (b) Corresponding plot for c-axis thermal conductivity. Counterintuitively, the correct 
use of the ab-axis MFPs in Eq. 6 leads to a higher thermal conductivity than using the c-axis 
MFPs, even though the ab-axis MFPs are larger than those along the c-axis. 

isotropic one. As a means of comparison, we instead calculate the thermal conductivity 
reduction using the longer ab-axis MFPs in the boundary scattering function rather than 
c-axis MFPs in the same hgure. The result shows the trend expected of an isotropic solid, 
with the thermal conductivity approaching zero for very thin dims and slowly approaching 
the bulk value as the him thickness increases. This prediction differs greatly from the correct 
calculation because of the form of Eq. 6: the thermal conductivity reduction depends only on 
the component of the MFP along the thickness direction, regardless of the ab-axis MFPs. 
As a result, even extremely thin graphite Elms maintain their high thermal conductivity 
despite the long MFPs along the ab-axis. 

The same calculation as above but for a thermal gradient along the c-axis is shown in Fig. 
2b. In this case, the graphite is inhnite along the c-axis and has a hnite thickness along one 
basal axis, similar to the geometry used in a prior experimental study.The thermal gradient 
exists along the c-axis of the graphite. The thermal conductivity versus him thickness is 
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FIG. 3: Spectral c-axis thermal conductivity versus phonon frequency using Eq. 6 (crosses) and 
with the c-axis MFPs in the boundary scattering function (circles). Even though the ab-axis MEPs 
are longer than those in the c-axis, the correct boundary scattering function only minimally affects 
phonons with long c-axis MFPs, resulting in a higher thermal conductivity than if the shorter 
c-axis MFPs are used in the boundary scattering function. 


presented according to Eq. 6, which uses the ab-axis MFPs, and for comparison the same 
calculation with the c-axis MFPs in the boundary scattering function. The calculation shows 
that, counterintuitively, the c-axis thermal conductivity obtained with the correct ab-axis 
MFPs in the boundary scattering function is actually higher than that with the c-axis MFPs, 
even though the ab-axis MFPs are considerably longer than those in the c-axis. The reason 
is that modes that primarily contribute to c-axis thermal conductivity have short in-plane 
MFPs and thus are minimally affected by the boundary function in Eq. 6. In contrast, 
using Eq. 6 with the c-axis MFPs incorrectly reduces the contribution of the long MFP 
c-axis modes that contribute the most to c-axis thermal conductivity, resulting in a smaller 
thermal conductivity than predicted by Eq. 6. 
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This point is illustrated in Fig. 3, which shows the spectral c-axis thermal conductivity 
^c,k = versus phonon frequency calculated using the ab- and c-axis components 

of the MFPs in the boundary scattering function S{^y). The correct boundary scattering 
function based on the ab-axis MFPs only minimally affects the phonons that primarily con¬ 
tribute to c-axis heat conduction because these phonons have short ab-axis MFPs and long 
c-axis MFPs. The primary effect of the correct boundary scattering function is to make the 
MFPs of phonons with small c-axis MFPs even smaller, thereby resulting in only a small 
reduction to the thermal conductivity. On the other hand, the incorrect boundary scatter¬ 
ing function signihcantly reduces the MFPs of phonons with long c-axis MFPs, resulting 
in a anomalously large reduction in thermal conductivity. These calculations show the im¬ 
portance of considering the dispersion anisotropies of a crystal when considering boundary 
scattering in thin dims. 

In summary, we have derived a generalized Fuchs-Sondheimer theory for crystals with 
arbitrary anisotropies and applied it to thin graphite films. We find that the large velocity 
anisotropy in graphite causes boundary scattering to have only a minimal effect on ab-axis 
thermal conductivity in films as thin as 10 nm. This observation has important implications 
for heat-spreading in electronic devices using thin graphite films and heat conduction in 
graphite foams. 

The author thanks Lucas Lindsay for providing the graphite dispersion and for useful 
discussions. This work was sponsored in part by the National Science Foundation under 
Grant no. CBET CAREER 1254213 and by Boeing under the Boeing-Caltech Strategic 
Research & Development Relationship Agreement. 
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